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Interpolation  from  Lagrange  to  Holberg 

Michel  Leger 


Abstract.  As  the  order  2 n tends  to  infinity,  Lagrange  interpolators  of 
periodically  sampled  ID  functions  converge  to  the  sine  function  modulated 
by  two  exponentials.  One  is  related  to  instabilities  and  the  other  to  Gaus- 
sian apodizing.  The  Hermite  interpolation  of  Lagrange  interpolators  gives 
convolutive  Cfc+1-differentiable  Lagrange-Hermite  interpolators.  Whereas 
their  support  has  width  of  order  2n  + 2,  the  active  part  of  their  impulse 
response  is  width  of  order  \/2n,  instead  of  2n  for  Holberg  interpolators, 
which  are  optimal  combinations  of  Lagrange-Hermite  interpolators,  and 
therefore  much  more  efficient.  Efficient  filters  can  be  derived  from  these 
differentiable  interpolators,  as  well  as  numerical  schemes  of  derivatives  at 
any  abscissa. 


§1.  Introduction 

Some  applications  require  very  large  2D  or  3D  regular  grids,  such  as  finite- 
difference  modeling  of  seismic  waves,  for  instance.  The  processing  of  these 
grids  involves  the  computing  of  numerical  schemes  of  first  or  second  deriva- 
tives, and  also  interpolators  and  filters.  These  quantities  need  to  be  evaluated 
very  efficiently  because  the  requirements  in  terms  of  computation  time  and 
memory  use  are  critical. 

Numerical  schemes,  interpolators  and  filters  are  interrelated  issues,  and 
I choose  to  study  them  from  the  viewpoint  of  interpolation,  which  is  the 
most  general.  For  sake  of  simplicity,  I assume  ID  interpolation  of  periodically 
sampled  functions,  with  unit-sampling  rate  and  even  orders. 

§2.  Lagrange  Interpolators 

Definition  1.  For  some  function  f with  known  values  fi  at  abscissae  Xi  = i, 
i € [— n,  n],  the  Lagrange  interpolator  ([9, 6, 2, 3])  of  order  2 n is 

4»({/*K=- »,*)  = £ f*  II  7TJ  = £ ^ P2n»-  (1) 

i——n  j — -n  J i=—n 
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It  can  be  shifted  and  centered  at  any  integral  abscissa  i: 

n 

4n({A}i'5-n.a')=  E 


3=-n 


Proposition  2.  As  2 n — > oo,  Lagrange  polynomials  converge  to  the  sine 
function  modulated  by  two  exponentials: 

:2 


Pt,i(X) 


5 «*(*)  exp(^-r)  (1  + 0{»-»)).  (2) 


r(z  — i) 


n + 


n + 2 


Proof:  Simple  changes  in  (1)  give 


(— l)2n— (n!)2a:  fV  _ ** 

P2n,i(x)  (n  -f  j)!  (n  - f)!  (x  - i)  ^ j 2 


(3) 


Since  n°°=i(l-f*) 


sin(7rrc) 


have 


ft  .(>-?)• 

j=l  J J=n+1 


(4) 


Since  In  n~  n+1(l~  fr)  = E“+i(f  + 0(J~4))>  and  since  J 
Ij-i/2  we  obtain 

II  (l-^)=exP(-^  + 0(n-3)). 


£ = + 


j=n+l 


(5) 


Moreover,  by  using  Stirling  formula  n!  — nn  exp(  n)  \/2im  (1  + 12„  + 288^  4" 
0(n-3)),  and  noting  that  ^(l+j^  + jgg^J  + OCn-'3))  = T2H+0(n  ) and  that 

ife  (2  - ^ - ^i)  = 0(n-3),  we  obtain  In  (n+^rr]T  = “(”  + *+  f)  ln(1  + 
i.)  _ (n  _ i + I)  ln(l  - i)  + 0(n-3).  Since  ln(l  +®)  = x-  ^-  + 3-  + 0(a:4), 
we  have  In  - 


(n+l)"fc^T  = + °(n'^  and  1161106 


(n!)2 


exp( — — y+0{n  3)). 


( n + i)!  (n  - i)\ 


n + 


(6) 


Noting  that  (-l^sinTraj  = sin7r(a:  - i)  and  inserting  (4),  (5)  and  (6)  in  (3) 
concludes  the  proof.  □ 

According  to  the  first  term  of  (2),  Lagrange  polynomials  converge  to  the 
perfect  interpolator  sine  function  (Fig.  1).  Away  from  the  center  of  the  interval 
of  the  data  points,  the  first  exponential  explains  the  well-known  instabilities 
(a  small  change  in  the  data  results  in  a large  change  in  the  interpolation),  and 
the  second  one  corresponds  to  Gaussian  apodizing  (see  apodization  in  [10]), 
that  is,  the  vanishing  effect  of  the  non-centered  data  points.  Note  that  these 
two  exponentials  compensate  one  another  as  x — ► i. 
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Fig.  1.  Lagrange  polynomials  of  order  14  divided  by  the  two  exponentials  of  (2), 
for  i = 0, 2, 4, 6.  Related  sine  functions  in  dotted  lines. 


§3.  Stationary  Lagrange  Interpolators 

Since  a Lagrange  interpolator  is  stable  only  near  the  middle  of  the  interval  of 
the  contributing  data  points,  a natural  idea  is  to  change  it  for  each  interval 
between  two  successive  points,  in  such  a way  that  the  interpolator  is  always 
used  at  its  best. 

Definition  3.  The  stationary  Lagrange  interpolator  of  order  2 n uses  at  ab- 
scissa x the  Lagrange  interpolator  of  order  2 n centered  at  ix , with  x = ix+dx, 

ix  € N,  4 € [0,  1[: 


^2n({/fc}fc£2Z,  x)  — ^2,  fix+i  P2n,j(dx)  — ^'2n(x)-  CO 

j——n 

Remark  4.  Any  interpolator  such  that  l({fk}k£2Z,x)  = Y2j=a  fix+i  Pj(dx), 
with  a,  b and  pj  independent  of  ix,  is  convolutive,  that  is,  there  exists  a 
continuous  function  A(x)  such  that  T({fk}kem,  x)  = A(f)*\(x),  with  A (/)  = 
/•  fi(x  ~ *)  being  the  “Dirac  comb”  modulated  by  function  f. 

This  is  obvious  by  considering  the  impulse  response  A(x)  = T({4o}fceZi  $)  = 
P-ix{dx),  with  Sk o the  Kronecker’s  symbol. 

As  a particular  consequence,  stationary  Lagrange  interpolators  are  con- 
volutive with  impulse  responses  A2 n(x)  — _ix  (dx).  Moreover,  they  are 
also  subject  to  Gaussian  apodizing  since,  as  2n  — > 00, 

, . . sin7rx  . —x2 

Mn  (x) exp  ( — — j- ) . (8) 

7 tx  n + 5 

It  is  clear  from  Fig.  2 that  these  impulse  responses  vanish  very  rapidly  as 
compared  to  the  sine  function. 

Remark  5.  Since  Vfc,  k € [1,  2n],  xk \2n(x)dx  = 0,  then,  with  T denot- 

ing the  Fourier  transform,  Vfc,  k £ [1,  2 n],  (^7(A2n)(w))„^o  = 0 ([!])- 

Hence,  stationary  Lagrange  interpolators  are  good  at  low  frequencies. 
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Fig.  2.  Impulse  responses  of  stationary  Lagrange  interpolators  for  orders  2 to 
14  (the  last  one  is  truncated),  sine  function  in  dotted  lines. 


§4.  Hermite  Interpolators 

In  a general  way,  Hermite  interpolators  are  consistent  with  the  values  of  a 
function  and  its  k (or  1 < i < l)  first  derivatives  at  l distinct  points  [4,2,3]. 
Here,  I consider  convolutive  interpolators  for  any  k and  l = 2. 

Definition  6.  The  Ck  piecewise  polynomial  Hermite  interpolator  is 

Tik{})  = A (/)  *Vko  + ■■■  + A(/W)  * r)kj  + . . . + A(/(fc))  * 1 7fcfc, 

where  r]kj(x)  are  [—1,  +l]-supported  basis  functions  such  that.  Vra,  m € [0,  fc], 
= 0 = Vk^\+l),  and  d?^0)  = Smj- 
From  now  on,  I consider  r]ko(x)  functions  only. 

Proposition  7.  The  polynomial  expression  of  Hermite  interpolators  r]ko  is 

j=2k  /i  y-k/nrj-* 

1 - mix)  = (2k  + 1)  ck2k  J2  1 -7— r xi+1  (9) 

U j + 1 

forx  € [0,  1],  andrjko(x)  = r)ko(—x)  Fora;  e [—1,  0].  Moreover,  asymptotically, 


[k 

lim77ko(a;)  ~ FlJ(x)  * 2y  — exp(— 4kx2),  when  k — > 00, 
with  n„(x)  = 1 if  x e [a  — |,  a + |]  and  n„(x)  = 0 otherwise. 


(10) 


Proof:  For  a;  € [0,  1],  it  is  clear  from  Def.  6 that  (1  — Tjko(x))^  = axk(l—  x)k . 
Since  7?ko(0)  = 1 and  r]ko(l)  = 0,  we  have 


- = [ (l-Vko{x))Wdx. 
J 0 


a 


(ii) 
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Fig.  3.  Right  part  of  impulse  responses  of  Hermite  interpolators  according  to 
(9),  their  asymptotic  form  according  (10)  in  dotted  lines. 

Let  us  define  Ikj  = Jq  xk(l  — x)ldx.  Integrating  by  parts  gives  Ik,  i = 
4+1, i-i,  that  is.  h,k  = (5^r4fc,o  = C%khk,  o by  recurrence.  Since 
Iik,o  = 2fc+I>  we  obtain  a = (2k  + l)Ckk,  and  therefore,  from  (11), 

(1  - W*))(1)  = (2fc  + 1)  Ckk  xk(l  - x)k.  (12) 

The  binomial  theorem  and  a simple  integration  concludes  the  proof  of  (9). 
Moreover,  by  the  change  of  variable  x = | +u  in  (12),  we  have  -77^(5  +11)  = 
(2k  + 1)  Ckk  4~k  (1  — 4 u2)k.  For  any  abscissa  u\  > 0 we  may  define  uk  = ^ 

and  we  have  exp(— 4fcu^)  = exp(— 4u2)  and  (1  — 4 uk)k  = (1  — 4^f)k.  Since 
lim*)^00(l  — 4^)fc  = exp(— 4«2),  then  (1  — 4u2)fc  ~ exp(—4ku2)  as  k — + 00. 
Moreover,  by  using  Ckk  = jjjyr  and  Stirling  formula,  we  obtain,  as  k — + 

00,  (2k  + l)Ckk4~k  ~ 2^|,  and  then  -J7$(x)  ~ 2yJ\  exp(-4fc(x  - §)2). 

Symmetrically,  we  have  ??^(x)  ~ 2yJ~^  exp(—4k(x  + |)2)  for  x € [—1,0]. 

Therefore,  ^^(x)  ~ exp(— 4 kx2)*  (S(x  + 1)  — S(x  — |)),  which  concludes 

the  proof  of  (10)  since  nj^(x)  = 5(x  + 5)  — <5(x  — ^).  □ 

Note  that  Lagrange  and  Hermite  interpolators  could  be  considered  as 
Fourier  pairs  since  the  right  members  of  (8)  and  (10)  are  mutual  Fourier 
transforms  if  7r2n  = 4k. 

§5.  Lagrange-Hermite  Interpolators 

Since  stationary  Lagrange  interpolators  are  good  at  low  frequencies  and  since 
Hermite  interpolators  are  differentiable,  hence  good  at  high  frequencies,  com- 
bining them  is  a natural  idea. 

Definition  8.  For  k > 0,  the  Lagrange-Hermite  interpolator  of  order  2 n is  the 
smooth  Hermite  interpolation  of  two  successive  Lagrange  interpolators,  that 
is, 


Mk+ 1(x)  = r)k0(x-ix)  Cfn(x)  + (1  - T)ko(x-ix))  £&+1(x). 
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Fig.  4.  Right  part  of  the  impulse  responses  of  (^-differentiable  Lagrange-Hermite 
interpolators  (a),  sine  in  dotted  lines.  Their  left  part  divided  by  the  sine 

_ 2 

(b),  as  compared  to  Gaussian  curves  exp(  in  dotted  lines. 

Proposition  9.  For  k > 0,  Lagrange-Hermite  interpolators  Mk„  1 are  Ck+1- 
differentiable. 

Proof:  Since  it  is  piecewise  polynomial,  1 's  C°°-differentiable  for  all 
non- integral  abscissas  x.  Moreover,  >s  convolutive,  and  thus  it  is 

sufficient  to  examine  it  around  x = 0.  From  Def.  8,  for  x > 0,  we  have 
= 7?)c0^2n+(1-%o)^2n-  Fro™  (12)>  we  have  7]k0{x)  = l-0(i*:+1),  and 
then  M^ix)  = £§„(a0  + 0(xk+1)  ( £\n  - £°n)-  Since  the  interpolators  £%n 
and  £\n  are  continuous,  £\n  - C%n  = O(x)  and  thus  1 = + 0(xk+2). 

Therefore,  for  any  j € [1,  k + 1],  A/f2n1^^(0+)  = £2iJ'(0+)-  The  same  argu- 
ment with  Mki 1 = Vko  £-2. n + (!  -Vko)  £%n  for  x < 0 leads  to  M ^ = 

£?2^\0~),  which  concludes  the  proof  since  £2„{x)  is  a polynomial.  □ 

The  Fourier  transform  of  Lagrange-Hermite  impulse  responses  (which  I 
call  pkn  for  interpolator  Mkn)  are  very  similar  to  those  of  stationary  Lagrange 
interpolators  below  the  sampling  frequency.  Beyond  the  sampling  frequency, 
for  the  orders  2 to  14,  the  rejection  in  dB  of  the  greatest  secondary  lobe  is 
27,  30,  32,  33,  34,  35  and  35  in  the  C°  case,  42,  47,  51,  53,  55,  57  and  58  in 
the  C1  case,  33,  36,  37,  38,  39,  40  and  41  in  the  C2  case,  and  slowly  decreases 
for  higher  differentiabilities.  From  this  viewpoint,  the  C1  Lagrange-Hermite 
interpolators  are  the  best  choice. 

§6.  Holberg  Interpolators 

Lagrange-Hermite  interpolators  are  good  at  low  and  high  frequencies,  but 
unsatisfactory  inbetween.  Indeed,  Gaussian  apodizing  makes  A2 n as  well  as 
pkn  gradually  ineffective,  since  their  cost,  which  is  proportional  to  the  length 
of  their  support,  increases  like  n,  whereas  their  active  part  widens  like  \/n. 
Faced  with  the  same  problem  in  terms  of  numerical  schemes,  Holberg  ([5])  had 
the  idea  of  combining  several  orders  and  optimizing  the  passband  for  given 
tolerance  and  maximal  order.  Holberg  interpolators  just  proceed  from  the 
same  idea. 
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Fig.  5.  Impulse  responses  of  (^-differentiable  and  0.1%-precise  Holberg  interpo- 
lators (a)  (Lagrange-Hermite  in  grey,  sine  in  dotted  lines).  Their  Fourier 
transform  (b)  and  an  enlargement  of  them  (c).  The  passband-ns-order 
diagram  (d)  showing  the  linear  behaviour  of  Holberg  interpolators  and 
the  parabolic  behaviour  of  the  Lagrange-Hermite  interpolators. 

Definition  10.  The  2n-hybrid  order,  Ck -differentiable  and  e-precise  Hol- 
berg interpolator  is  the  linear  combination  = YllZi  Pt^2i  such  that 

" Pi  = 1 and  such  that  the  following  passband  is  maximized, 

B(Pt,...,Pen-1)  = sup({«/;Vf<V;|^/af^)(0-l|  < *})• 

i 

Clearly,  these  Holberg  interpolators  are  convolutive,  with  impulse  responses 

thn  = El  Pitin- 

Fig.  5 illustrates  the  case  of  e = 0.1%  and  k = 1 (C1-differentiability). 
Fig.  5a  shows  that  the  impulse  responses  of  Holberg  interpolators  decrease 
much  slower  than  the  Lagrange-Hermite  interpolators.  The  sum  of  the  Pf  is 
one,  but  their  absolute  sum  may  be  far  from  one,  for  instance  about  4300  in 
the  case  of  the  14th  order.  In  the  Fourier  domain,  Fig.  5b,  together  with  its 
enlargement  Fig.  5c,  shows  that  the  passband  is  increased  a lot  from  Lagrange- 
Hermite  to  Holberg. 

For  tolerances  1%  and  0.1%,  Fig.  5d  shows  the  parabolic  behaviour  around 
infinite  order  of  the  passbands  of  the  Lagrange-Hermite  interpolators  (grey 
arrows  are  parabolas  with  vertical  tangent  at  the  corner).  This  is  due  to 
Gaussian  apodizing,  since  the  Fourier  transform  of  (8)  gives  the  convolution 
of  a box  function  with  a Gaussian  function  that  narrows  like  -4=  as  n — ► oo. 
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Fig.  6.  Right  part  of  the  impulse  responses  of  the  lOth-order  Lagrange-Hermite 
(in  grey)  and  Holberg  interpolators  of  first  (a)  and  second  (b)  deriva- 
tives. Perfect  (but  truncated)  interpolators  for  derivatives  are  displayed 
in  dotted  lines. 

On  the  contrary,  the  passbands  of  the  Holberg  interpolators  have  a linear 
behaviour  (black  arrows  are  straight),  which  is  the  best  possible  one  because 
the  active  part  of  the  impulse  response  of  the  interpolator  cannot  expand  faster 
than  its  support.  The  passbands  are  referred  to  the  sampling  frequency. 

The  optimization  of  passband  B in  the  (3f  space  is  difficult  because  B 
is  discontinuous  along  n — 1 hypersurfaces  (related  to  tangencies  at  lie) 
which  intersect  at  the  optimum.  To  overcome  this  difficulty,  I used  a method 
consisting  of  the  following  steps: 

1)  choose  a strictly  increasing  sequence  of  frequencies  i/;  in  ]0,  0.5 [,  with 
1 < i < n — 1, 

2)  set  at  1 + e the  value  of  the  combination  at  frequency  1 — s at 
1 + e at  i/n_3,  and  so  on  until  iq,  and  finally  1 at  frequency  j'o  = 0, 

3)  solve  the  linear  system  for  the  0f, 

4)  detect  the  frequencies  at  which  the  combination  is  extremal, 

5)  if  these  frequencies  are  close  to  the  i then  stop,  else  update  the  and 
come  back  to  step  2. 

A priori,  the  feasability  of  steps  3 and  4 and  the  convergence  are  not 
guaranteed.  In  practice  however,  this  method  works  simply  well,  with  less 
than  ten  iterations.  See  also  [7,8]  on  optimal  filtering. 

§7.  Applications 

The  main  application  of  Holberg  interpolators  are  Holberg  numerical  schemes 
([5])  because  they  are  cost-effective  in  the  field  of  numerical  simulation  of 
acoustic  wave  propagation.  Especially  in  3D,  this  effectiveness  is  of  consider- 
able importance  because  of  the  huge  amount  of  computing  time  needed. 

From  a 2n-order  C^-differentiable  Lagrange-Hermite  or  Holberg  interpola- 
tor, 2n-order  Cfc_m-differentiable  interpolators  of  the  mth  derivative  (m  < k ) 
can  be  easily  derived,  as  well  as  numerical  schemes  of  these  derivatives  at  any 
abscissa.  Fig.  6 shows  the  responses  of  the  lOth-order  Lagrange-Hermite  and 
Holberg  (e  = 0.1%)  interpolators  of  first  (a)  and  second  (b)  derivatives.  The 
values  at  integral  abscissas  give  standard  numerical  schemes.  In  a similar  way, 
the  integration  of  these  impulse  responses  could  result  in  Newton-Cotes-like 
Holberg  formulas  (see  Newton-Cotes  Formulas  in  [10]). 
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0 abscissa  (km)  2.775 


Fig.  7.  A seismic  section  (a)  has  been  horizontally  filtered  (b)  for  antialiasing, 
threefold  undersampled  (c),  and  finally  interpolated  (d). 

Holberg  interpolators,  and  numerical  filters  that  can  be  generated  from 
them,  are  also  interesting  for  their  efficiency.  For  instance,  the  response  of 
a 1%-precise,  6th-order,  C1  Holberg  interpolator  has  been  fourfold  oversam- 
pled and  used  to  filter  horizontally  Fig.  7a  into  Fig.  7b  for  antialiasing.  The 
threefold  undersampling  of  Fig.  7b  gives  Fig.  7c.  The  1%-precise,  6th-order, 
C1  threefold  Holberg  interpolation  of  Fig.  7c  gives  Fig.  7d.  The  similarity  of 
(b)  and  (d)  measures  the  quality  of  the  filtering  and  of  the  interpolation. 

§8.  Conclusions 

In  the  case  of  periodic  data  points,  Lagrange  interpolators  converge  to  a sine 
function  multiplied  by  two  exponentials.  The  first  one  explains  the  well- 
known  instabilities  of  Lagrange  polynomials,  which  only  vanish  at  the  center 
of  the  interval  of  the  data  points.  The  second  exponential  explains  the  van- 
ishing influence  of  non-centered  data  points  (Gaussian  apodizing).  Station- 
ary Lagrange  interpolators  are  stable  and  convolutive.  Hermite  interpolators 
are  as  differentiable  as  desired  and  convolutive  and  their  impulse  response 
converges  to  the  convolution  of  a box  function  with  a Gaussian  function. 
Lagrange-Hermite  interpolators  combine  the  advantages  of  unlimited  order 
and  differentiability.  Because  of  Gaussian  apodizing,  these  interpolators  be- 
come ineffective  at  high  orders.  On  the  other  hand,  Holberg  interpolators 
have  a much  better  quality/cost  ratio  since  they  are  optimal  combinations  of 
Lagrange-Hermite  interpolators.  From  Holberg  interpolators,  efficient  numer- 
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ical  schemes  of  derivatives  can  be  evaluated  at  any  abscissa,  and  oversampling 
their  impulse  response  gives  short  but  efficient  filters. 
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